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ABSTRACT 



Realistic 3D radiative MHD simulations reveal the magneto- convective processes underlying the formation of the pho- 
tospheric fine structure of sunspot s, including penumbral filaments and umbral dots. Here we provide results from 
a statistical analysis of simulated umbral dots and compare them with reports from high-resolution observations. A 
multi-level segmentation and tracking algorithm has been used to isolate the bright structures in synthetic bolometric 
and continuum brightness images. Areas, brightness, and lifetimes of the resulting set of umbral dots are found to be 
correlated: larger umbral dots tend to be brighter and live longer. The magnetic field strength and velocity structure of 
umbral dots on surfaces of constant optical depth in the continuum at 630 nm indicate that the strong field reduction 
and high velocities in the upper parts of the upflow plumes underlying umbral dots are largely hidden from spectro- 
polarimetric observations. The properties of the simulated umbral dots are generally consistent with the results of recent 
high-resolution observations. However, the observed population of small, short-lived umbral dots is not reproduced by 
the simulations, possibly owing to insufficient spatial resolution. 
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1. Introduction 

Recent realistic MHD simulations have revealed that the 
complex brightness and flow structure of sunspots re- 
sults from convective energy transport dominated by a 
strong vertical (in the umbra) or inclined (in the penum- 
bra) magnetic field ( [Heinema nn et al.ll2007l : IScharmer et al.1 
2008; R empel et al.| [2009b,a). In the umbra, such magneto- 
convection occurs in t he form of narrow upflow plumes with 
adjacent downflows (iSchiissler fc Voglerl 120061 henceforth 
referred to as SV06). In computed intensity images, these 
plumes appear as bright features that share various proper- 
ties with observed umbral dots. Here we provide a system- 
atic study of the umbral simulation results with specific 
emphasis on the properties of the bright structures and 
their comparison with recent results from high-resolution 
observations. 

2. MHD simulation 



We u sed the MURaM code ([Voderl 120031 : IVogler et all 
l2QQ5h with nearly the same simulation setup as SV06. 
The dimensions of the computational box are 5.76 Mm x 
5.76 Mm in the horizontal directions and 1.6 Mm in the ver- 
tical. Rosseland optical depth unity is about 0.4 Mm below 
the upper boundary. The grid cells have a size of 20 km 
in the horizontal directions and 10 km in the vertical. The 
magnetic diffusivity is 2.8-10 6 m 2 s _1 , the lowest value com- 
patible with the grid resolution. The hyperdi ffusivities for 
the o ther diffusive processes (for details, see IVogler et al.l 
2005) in the deepest layers of the box were chosen lower 
than the values used by SV06 in order to minimize their 
contribution to the energy transport. Since we do not syn- 
thesize spectral lines or Stokes profiles, the radiative trans- 



fer can be treated in the gray approximation. The fixed 
vertical magnetic flux through the computational box cor- 
responds to a mean field strength of 2500 G. Side bound- 
aries are periodic, and the top boundary is closed for the 
flow. The lower boundary is open and the thermal energy 
density of the inflowing matter has been fixed to a value of 
3.5 • 10 12 erg-cm -3 , leading to a total surface energy output 
comparable to that of a typical sunspot umbra, i.e., about 
20% of the value for undisturbed solar surface. The mag- 
netic field is assumed to be vertical at the top and bottom 
boundaries. 

We started our simulation from the last snapshot of the 
run analyzed by SV06 and continued for 10.8 hours solar 
time in order to obtain a sufficiently large statistical sample 
of simulated umbral dots. Since the average magnetic field 
in the simulation is vertical, our results represent the inner 
part of a sunspot umbra with central umbral dots - as op- 
posed to peripheral umbral dot s, which often are related t o 
penumbral filaments (e.g.. iGrossmann- Doerth et aflll986D . 



3. Image segmentation 

To obtain statistical properties of the large number of sim- 
ulated umbral dots (UDs), we used an automated detec- 
tion method to distinguish between bright features and 
the surrounding dark background (DB). Since both the 
UDs and the background cover a broad range of in- 
tensities, we chose the MLT (Mult i Level Tracking) al- 
gorithm of iBovelet fc Wiehrl (|200ll ) as the basis of our 
image segmentation. This method has been successfully 
appli e d to various observatio n al da t a (e.g.. IWiehr et all 
120041 : IBovelet fc Wiehrl 120071 . 120081 : iKobel et al.1 12009). 
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Fig. 1. Over-segmentation of simulated umbral dots due to 
dark lanes. The section of the original image (middle panel) 
shows umbral dots with dark lanes. The unmodified MLT 
algorithm leads to an undesired splitting of many dots in 
the segmented image (left panel). Introducing a modifica- 
tion of the algorithm controlling the merging of segments 
removes the problem in most cases (right panel). 




Fig. 2. Result of the MLT image segmentation. Left: orig- 
inal image; Middle: segmented image (binary map); Right: 
segmented umbral dot areas with their original intensities. 



i ncluding high-resolution images of sunspot umbrae 
(jRiethmuller et alJl2QQ8bh . 

The basic concept of the MLT algorithm is to segment 
the image by applying a sequence of intensity thresholds 
with decreasing values, keeping all features already de- 
tected. In our case, a complication arises from the dark 
lanes shown by many of our simulated UDs. This can lead 
the algorithm to split an UD into two parts (see Fig. [1]). 
In order to avoid such undesired splitting, we permit the 
merging of segments that were separate for threshold n after 
applying threshold level n+l provided that their maximum 
intensity does not exceed a value equal to the latter thresh- 
old multiplied by a suitably chosen factor c. For our analy- 
sis, we used 15 intensity thresholds ranging from about 2.5 
to 1.5 times the intensity of the dark background (some- 
what depending on the dataset used, see Sec. 4.2) and a 
value of c = 1.45, which suppresses most cases of unwanted 
splitting of UDs. 

In a final step of the segmentation process, we removed 
from the list of UDs all segmented structures whose max- 
imum intensities do not exceed the lowest threshold by at 
least 10%; such features are considered to be mere fluctu- 
ations of the background. Fig. [2] shows an example of the 
resulting segmentation of a typical simulated image. 

In order to study UD evolution and lifetimes, we ex- 
tended the MLT algorithm to include time, such that the 
segmentation can be carried out in three dimensions (two 
spatial and one temporal) to follow the temporal evolution 
of individual UDs. 



4. Results 

4.1. Overall temporal evolution 

The relative inefficiency of the convective energy transport 
in the strong umbral magnetic field leads to a rather long 
thermal relaxation time of the system, so that it can take 
hours until it settles down into a statistically stationary 
state. With our long simulation run of 10.8 hours solar time, 
we were able to follow the overall changes of the system 
throughout the relaxation phase. 

The steady decrease of the total internal energy in the 
computational box shown in Fig. [3] (left panel) clearly in- 
dicates that the sytem is not relaxed until about 6-7 hours 
after the start of the simulation. Comparing the time- 
integrated total radiative output from the box with the 
integrated decrease of the total internal energy given in the 
right panel of Fig. [3] suggests that during the first half of the 
evolution the energy output was mostly covered by a loss 
of internal energy, i.e., the plasma in the box cooled. Only 
after about 6 hours, the two curves start to diverge and the 
internal energy does no longer decrease. At about the same 
time, the slope of the integrated energy output steepens, 
indicating a somewhat higher energy output. This is con- 
firmed by Fig. 2] (left panel), which shows that the total 
radiative flux is higher by about 10% in the second half of 
the run. 

These results suggest that the character of the energy 
supply has changed halfway through the run: while the ra- 
diated energy was taken from the internal energy in the 
first half and was not significantly replenished through con- 
vective inflows through the bottom of the box, convection 
became more effective in the second half of the run and led 
to a thermally relaxed stationary situation with a some- 
what higher total radiative output. This interpretation is 
supported by considering the time evolution of kinetic en- 
ergy (Fig. 2J right panel), which shows an increase by a 
factor of 3 in the second half of the run. The magnetic en- 
ergy also grows slightly, by about 10%, probably due to the 
stronger fluctuations (local compressions of the magnetic 
flux) caused by the higher velocities. 

4.2. Statistical properties of umbral dots 

The existence of the two phases (cooling phase and con- 
vective phase) in the temporal evolution of the umbra sim- 
ulation is also reflected in differences of the average UD 
properties in both phases. The snapshots of (bolometric) 
brightness in Fig. [5] indicate that larger UDs appear in the 
second phase while the total number of UDs decreases. We 
therefore defined separate datasets for each phase and ana- 
lyzed them individually. The first part (PI) covers the time 
between 0.5 h and 4.5 h after the beginning of the run while 
the second part (P2) extends from 6.2 h to 10.3 h. We ex- 
cluded the transition between the two phases and only kept 
periods for which the average UD properties do not show 
significant secular trends. 

For the MLT segmentation considered in this section, we 
used maps of the bolometric (frequency-integrated) emer- 
gent intensity in the vertical direction. The maps are sep- 
arated in time by about 1 min each, so that UDs can be 
analyzed in a statistical manner as well as their evolution 
followed in time. Owing to the different characteristics of 
the UDs, slightly different MLT threshold levels were chosen 
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time [h] time [h] 

Fig. 3. Thermal relaxation of the system. Left panel: time evolution of the total internal energy (normalized to its initial 
value); Right panel: total (radiative) energy output (solid line) and loss of total internal energy (dashed line), both 
integrated in time from t = onward. 




time h 




time h 



Fig. 4. Temporal evolution of the system. Left panel: radiative energy flux from the box, normalized by the corresponding 
value for the same area of undisturbed solar photosphere; Right panel: total kinetic energy (solid line) and magnetic energy 
(dashed line), both normalized to their values at t = 0. 



for PI and P2, respectively. For dataset PI, the 15 inten- 
sity thresholds range between 20% and 31% of the average 
quiet-Sun brightness while for P2 they range from 22% to 
35%. 

Fig. [6] shows, separately for PI and P2, the temporal 
evolution of various average quantities of the segmented 
UDs and of the remaining domain area, the DB. The area 
fraction of UDs (top left panel) tends to be somewhat 
higher in P2 compared to PI; its relatively strong fluc- 
tuations result from the relatively small size of the com- 
putational box. In P2, there are less (top right panel) but, 
on average, larger (middle right panel) and brighter (middle 
left and bottom left panel) UDs, as already indicated by the 
snapshots shown in Fig. [5l The DB is also slightly brighter 
in P2 (bottom left panel). The mean bolometric intensity of 
the UDs approaches roughly stationary values of about 0.23 
(PI) and 0.26 (P2) of the quiet-Sun value, or about 1.6 (PI) 
and 1.7 (P2), respectively, in units of the corresponding DB 
intensity (bottom right panel). These results indicate that 
the upflows underlying the UDs contribute more strongly 




012345 012345 
[Mm] [Mm] 



Fig. 5. Snapshots of the bolometric emergent intensity 
from the two phases of the simulation. Upper panel: ini- 
tial cooling phase (t ~ 2h); Lower panel: convective phase 
(t~9h). 



to the convective energy transport in P2. This includes not 
only the fraction of the radiative output directly covered 
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Fig. 6. Time evolution of UD properties, given separately 
for the two phases PI (cooling phase, red curves) and P2 
(convective phase, green curves). A: UP area fraction; B: 
number of UDs; C: average UD area; D: mean bolometric 
UD intensity (normalized by the mean quiet-Sun value); 
E: average bolometric intensity of DB (normalized by the 
mean quiet-Sun value); F: mean bolometric UD intensity 
(normalized by the corresponding average for the DB). 
Time is counted from the start of the respective dataset 
(0.5 h and 6.2 h after the start of the simulation for PI and 
P2, respectively.) 



by the UDs (about 30%), but also the horizontal radiative 
losses of the upflow plumes, which heat their environment 
and thus contribute to the energy output of the DB. The 
higher average brightness of the DB in the second phase 
(P2) is thus consistent with a bigger contribution of the 
UDs to the overall convective energy transport. 

The average UD area resulting from our simula- 
tion and segmentation analysis is 0.08 Mm 2 for PI and 
0.14 Mm 2 for P2, corresponding to average diameters of 
320 km and 420 km under the assumption of circular 
areas. These values are higher than those reported in 
recent observa tional studies of UDs w ith the 1-m SST 
on La Palma: IRiethmull er et aD (l2008bT) give a n average 
area of 0.04 Mm 2 , ISobotka fc H anslmeier (2005) estimate 
0.025 Mm 2 while ISobotka & : Puschmand pOQ9) even find 
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Fig. 7. Area histograms of bolometric intensity (upper 
panel) and vertical field strength (lower panel). The field 
strength is given taken at z = 0, corresponding to the av- 
erage geometrical height level of Rosseland optical depth 
unity. Solid curves refer to UDs, dashed curves to DB. Red 
and green curves indicate datasets PI and P2, respectively. 



a value near 0.01 Mm 2 . On the other hand, the UDs with 
dark lanes studied by the latter authors have an average 
area a round 0.08 Mm 2 and even la rger UDs have be e n stud- 
ied bv lBharti etafl (|2007b[ ). Also. lWatanabe et all (|2009aD 
find an average area of about 0.1 Mm 2 from Hinode data. 
In such a comparison one has to take into account that the 
lowest intensity thresholds in our MLT segmentation are 
not far above the DB values, while the observational anal- 
yses typically define the edges of an UD by the half-width 
of the local intensity contrast or by the inflection point of 
the intensity profile. These procedures tend to yield smaller 
areas than our approach, which covers the full extension 
of the UDs (cf., for example, the UP o utlines shown in 
Figs. 5 and 6 of IRiethmuller et al.ll2008bl). We tested this 
conject ure by applying the procedure of IRiethmuller et al.l 
(|2008bf ) to a few representative snapshots in our datasets. 
We found that the resulting UP areas are, on average, 
about 50% smaller than the values determined by our MLT 
procedure. This brings the simulation and the observa- 
tional results into rough agreement - the simulations pos- 
sibly lacking a population of very small, short-lived UPs 
([Sobotka fc Puschmannll2009T ). 

For the full datasets PI and P2, Fig. shows area his- 
tograms of the bolometric intensity and the vertical compo- 
nent of the magnetic field at constant height z = 0, which 
corresponds to the mean level of Rosseland optical depth 
unity. Histograms for UPs and PB are given separately (in- 
cluding all image pixels belonging to each class). The UP 
intensity distributions reflect the fact that UPs tend to be 
brighter in P2. The tail of PB intensities overlapping with 
the UP distribution results from the 'trenches' generated 
by the MLT algorithm to separate closely neighboring or 
clustered UPs and the exclusion of the faintest UPs. 

The area histograms of the vertical magnetic field given 
in the lower panel of Fig. shows a very broad distribution 
of field strength in the UPs, ranging from slightly negative 
values up to about 3000 G. While the field strength becomes 
very low owing to expansion and flux expulsion in the core 
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Fig. 8. Histograms of UD area (left) and mean bolometric 
brightness (right, normalized to the average DB) for all 
segmented maps. Red and green curves refer to datasets 
PI and P2, respectively. 



of the rising plumes that generate the UDs, horizontal ra- 
diative heating extends the wings of the bright intensity 
structure (i.e., the UD) into the surrounding umbra with 
high field strength. Since we segment the images according 
to the intensity structure, the peripheries of the so-defined 
UDs harbor strong magnetic fields. Negative field strength 
arises when a strong downflow catches and drags a mag- 
netic field line, so that it is stretched into a hairpin-like 
shape (see Fig. [14] and Sec. 4.3). The average field strength 
in the DB is higher in P2 since the larger area fraction of 
UDs compresses the magnetic flux in the DB. 

The distributions of UD area and mean brightness for all 
segmented maps (with one minute cadence, so that UDs are 
considered more than once, in various stages of their evo- 
lution) are illustrated by the histograms shown in Fig. [HJ 
UDs with areas below 0.02 Mm 2 were omitted as most of 
them represent fluctuations in the DB. The area distribu- 
tion for P2 is broader, with less smaller and more larger 
UDs than for PI. The largest UDs have areas of 0.44 Mm 2 
and 0.79 Mm 2 for PI and P2, respectively. The distribution 
of average UD bolometric brightness contrast (right panel 
of Fig. [HJ) is shifted towards higher values for P2, the mean 
values being 1.62 (PI) and 1.72 (P2). 

Since the MLT tracking algorithm identifies and tracks 
UDs during their whole life cycle, we determined lifetimes 
for 558 UDs from PI and 500 UDs from P2. The life- 
time distributions are shown in the form of histograms 
in Fig. [9] (excluding the roughly 10% of the UDs with 
lifetimes below 8 minutes, which seem to constitute a 
separate group of low- amplitude fluctuations). The life- 
time distributions are similar for both datasets, with av- 
erage values of 28.2 min (PI) and 25.1 min (P2). The 
mean lifetimes reported from observations range from 
3 min to about an hour, with recent high-resolution 
observations typically indicating rather short mean life- 
times below 10 minutes (e.g., Kusoffskv fc Lu ndstedt 
1986 t lEwelll Il992| : ISobotka et ail Il997l: iKitai et all 12007: 
Riethmuller et all l2QQ8bl : ISobotka fc Puschmannl I2OO9I : 
Wa tanabe et al.ll2009al ). However, similar to our simulation 
results, the distribution of lifetimes from any observation is 
rather broad and 'typical' l ifetimes cannot be defined easily 
([ So cas- Navarro et~ al. 2004). Taken at face value, the obser- 
vations indicate the existence of a large number of short- 
lived UDs that are not shown by the simulations. However, 
it is also possible that some observational lifetimes are af- 
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Fig. 9. Histogram of lifetimes for UDs that could be 
tracked during their whole evolution. Red and green curves 
refer to datasets PI and P2, respectively. 
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Fig. 10. Relationship between UD lifetime, area and 
brightness. Left: scatter plot of lifetime vs. maximum area; 
Right: scatter plot of lifetime vs. maximum area- aver aged 
bolometric brightness reached during the lifetime of an UD 
(normalized by the average brightness of the DB). Red and 
green dots refers to datasets PI and P2, respectively; solid 
black (PI) and blue (P2) lines connect binned values for 50 
points each. 



fected by brightness fluctuations of UDs and by seeing ef- 
fects. 

The relationship between lifetime and maximum UD 
area and brightness in the course of the UD evolution is 
shown in Fig. [10] in the form of scatter plots with curves 
connecting binned values. Although there is a significant 
amount of scatter, the plots indicate that longer-lasting 
UDs tend to be larger and brighter. A qualitatively sim- 
ilar life time-size relation is reported by iRiethmuller et aD 
(2008b) for UDs with lifetimes below 20 min while longer- 
lived UDs are not found to be systematically larger. 

The lifetime-area and lifetime-brightness correlations 
are probably due to the fact that stronger and more ex- 
tended convective upflows are maintained longer and ere- 
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Fig. 11. Relation between area and brightness for the same 
set of UDs as in Fig. [10l Left: scatter plot of maximum 
area vs. maximum (normalized) bolometric brightness dur- 
ing the lifetime of an UD; Right: scatter plot of mean UD 
area (over UD lifetime) vs. mean area- averaged brightness 
(over UD lifetime). Red and green color refers to datasets 
PI and P2, respectively; solid black (PI) and blue (P2) 
lines connect binned values for 50 points each. 




Fig. 13. Profiles of vertically emerging intensity along a 
cut at x = 4.6 Mm through the image shown in the right 
panel of Fig.[5l The solid line refers to the continuum inten- 
sity at 630 nm, the dashed line to the bolometric intensity. 
Intensities are normalized individually by the correspond- 
ing local DB intensities. 
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Fig. 12. Area histograms of the continuum intensity at 630 
nm. Solid curves refer to UDs, dashed curves to DB. Red 
and green curves indicate datasets PI and P2, respectively. 



ate larger and brighter UDs: higher upflow speed entails 
a bigger convective energy flux density (brighter UDs) 
and also more mass flux and kinetic energy available for 
the sideways expansion of the upflow plume (bigger UDs). 
This explanation is in line with the relationship between 
UD area and (maximum and mean) brightness shown 
in Fig. ITT1 Similar t r ends were fo und observationally by 
iTritschler fe Schmidtl (l2002h an d bvlSobotka fc Puschmannl 
(|20Q9h while IRiethmuller et all j2008b|) do not find a clear 
relationship. 

4.3. UD properties and optical-depth dependence at 630 nm 

The iron line pair near 630 nm and the nearby continuum is 
often used for observations of sunspot fine structure. In or- 
der to compare with observational results derived from such 
data, we considered 57 (PI) and 74 (P2) snapshots from our 
simulation and calculated images of the (vertically emerg- 
ing) continuum intensity at 630 nm. These images were 
segmented by the MLT procedure to obtain maps of UDs. 
In addition to the UD properties studied in the previous 
section, we correlated the UD maps with the distributions 
of vertical magnetic field and vertical velocity on surfaces of 
equal optical depth at 630 nm. A full synthesis of line pro- 
files and Stokes parameters for a detailed comparison with 
observed spectro-polarimetric data is beyond the scope of 
this work and will be carried out in a forthcoming paper 
(Vitas, Vogler & Keller, in preparation). 



Area histograms for the continuum intensity at 630 nm 
are shown in Fig. [T2J Compared to the histograms of the 
bolometric intensity shown in Fig. [71 the intensity contrasts 
between UDs and DB is much higher, resulting from the 
fact that 630 nm is on the blueward side from the maxi- 
mum of the respective Planck functions. The significantly 
higher contrast between DB and UD is also illustrated by 
the intensity profiles shown in Fig. [13l The mean UD in- 
tensity values in Fig. [12] are 2.58 for PI and 2.88 for P2 
(both relative to the corresponding average DB). These val- 
ues are consistent with the i ntensity ratios reported from 
observations at 602 nm (e.g. , ISobotka fc Hanslmeierl l2005; 
ISobotka fc Puschmannl [2009). Individual UDs in the simu- 
lation can reach significantly higher intensities (see Fig.[T3|). 

Maps of the vertical magnetic field and vertical velocity 
on surfaces of equal optical depth at 630 nm are given in 
Fig.[H While the magnetic field distribution becomes more 
diffuse and homogeneous with height, the velocity remains 
rather intermittent. Strong downflows at the edges of the 
largest UDs at T630 = 1 (left panels) can capture inclined 
magnetic field lines and drag them downwards, thereby cre- 
ating a hairpin-like structure with a patch of reversed po- 
larity, indicated by red color in the magnet ic- field map. 
Quantitative information on the distributions is provided 
by the area histograms for the same optical-depth levels 
shown in Fig. [15l The histograms at T630 = 0.01 roughly 
represent the results that would be obtained by carrying out 
inversions of spectro-polarimetric observations in the neu- 
tral iron lines at 630.15 nm and 630.25 nm. The histograms 
at T630 = 1 are roughly similar to those at constant geomet- 
rical depth z = shown in the lower panel of Fig. [7J with 
a broad range of values in the UDs. At tqso = 0.1, field 
strengths below 1000 G are rarely found in UDs (particu- 
larly for those from PI) and at T630 = 0.01, the UD dis- 
tributions are further shifted to higher field strengths. For 
instance, 90% of the UDs of dataset PI have field strengths 
between 1800 G and 2200 G at t 630 = 0.01. This effect 
is due to the elevation of the surfaces of constant optical 
depth in UDs (in comparison to the DB) and the cusp-like 
shape of the UP structure (cf. I Schussler & Vogler 2006; 
IRiethmuller et all l2QQ8aL iBharti et all l2QQ9h . The larger, 
stronger UDs in the convective phase (dataset P2, green 
lines) reach higher into the umbral atmosphere [cf. Fig. [18] 
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Fig. 14. Maps of vertical magnetic field strength (upper row) and vertical velocity (lower row) on surfaces of constant 
optical depth at 630 nm (from left to right: T630 = 1, 0.1, 0.01, respectively). The color table for the magnetic field ranges 
from 3780 G (blue) over zero (white) to —190 G (red). For the velocity the range is from 1.9 kms -1 (blue, upflow) over 
zero (white) to —1.7 kms -1 (red, downflow). 



and show therefore still a stronger signature of the reduced 
field strength at T630 = 0.01 than the UDs from PI. For 
both datasets, the distributions for UDs and DB approach 
each other as the optical depth decreases; this reflects the 
decreasing area fraction of the cusp-shaped UD structures 
with increasing height. 

The corresponding area histograms of the vertical ve- 
locity given in Fig. [16] show rather weak motion in the 
DB: at all three r levels, about 50 % of DB shows al- 
most no vertical flow while the rest of the area has veloci- 
ties in the range ± 200 ms -1 . The strong upflows creating 
the UDs have already been considerably braked when they 
reach the level T630 = 1. Higher up at levels T630 = 0.1 
and T630 = 0.01, high- velocity tails (up to ~ ilkrns -1 ) 
appear in the distributions. These are mostly due to the 
upward- directed, jet- like Valve flows' from the UD cusps 
(|Choudhurilll986l : ISchiissler fc Voderll2006t ) and the corre- 
sponding return flows. The latter, however, may be affected 
by the presence of the closed upper boundary of the simu- 
lation box. 

Figure fT71 illustrates, by means of histograms, the prop- 
erties of average vertical magnetic field and velocity at var- 
ious optical depths, together with the average elevations of 
the optical-depth levels for all UDs from the datasets PI 
and P2. The magnetic field histograms show that, at all 
optical-depth levels, the expulsion of magnetic flux from 



the expanding upflow plumes underlying the UDs is more 
clearly reflected in the case of P2. However, in most cases 
the region of very low field strengths does not extend above 
the surface of tqso = 1- The cusp shape of the ve rtical UD 
structure (see Fig. 2 of Schuss ler fc Voglerl 120061 ) leads to 
an increase of the field strength for the iso-r surfaces above 
T630 = 1. In addition, the UD area is defined by the inten- 
sity structure corresponding to a height near T630 = 1; the 
cusp's cross section shrinks toward higher iso-r surfaces, so 
that we progressively sample more of the strong-field re- 
gion surrounding the plume. A similar situation is found 
for the upflows (second row), where the average velocities 
do not exceed a few hundred ms -1 at tqso = 1. At higher 
levels, the jet-like outflows along the cusps of the upflow 
plumes are seen. The distribution of upflows is shifted to- 
wards lower velocities for the UDs from P2, since bigger and 
flatter UDs show weaker outflows. This results from the fact 
that the height reached by the upflow plume is only weakly 
dependent on its size (since the stratification above optical 
depth unity is strongly subadiabatic), so that the cusp over- 
lying a large UD is broader (has a larger aspect angle) and 
the whole structures thus is flatter. Therefore, the chan- 
nelling of the upflowing matter into the top of the cusp is 
less efficient and the jet-like outflows are weaker. The down- 
flows (third row) at the UD periphery (for T630 = 1) and 
adjacent to the upflow jets (for the other levels) are similar 
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Fig. 15. Area histograms of the vertical magnetic field 
component on surfaces of constant optical depth at 630 nm. 
Upper panel: tq^q = 1; middle panel: tqsq = 0.1; lower panel: 
7630 = 0.01. Solid lines refer to UDs, dashed lines to the 
DB. The red and green colors indicate datasets PI and P2, 
respectively. 
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Fig. 16. Same as Fig.[l5]but for the vertical component of 
the velocity. Positive velocity corresponds to upward mo- 
tion, negative velocity to a downward flow. 
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Fig. 17. Distributions of averages over UD area of vertical 
magnetic field (top row) and vertical velocity (separately for 
upflows, v z > 0, and downflows, v z < in the second and 
third row, respectively), all on surfaces of constant optical 
depth at 630 nm (left column: tqso = 1, middle column: 
7"630 = 0.1, right column: tqso = 0.01). The bottom row 
gives the elevation of the optical-depth levels of the UDs 
with respect to the corresponding mean levels in the DB. 
Red and green lines refer to UDs from datasets PI and P2, 
respectively. 



in average magnitude for PI and P2 while the distributions 
for P2 are somewhat broader. 

The bottom row of Fig. [TT] shows histograms of the rela- 
tive height difference between the levels of constant optical 
depth in UDs (averaged over UD area) and in the mean DB, 
thus representing the elevation of the optical-depth levels 
in the UDs relative to the DB. While the mean elevations 
range between 60-90 km, the levels for the UDs from P2 
are, on average, 10-20 km higher than the corresponding 
levels of UDs from PI. Thus, the P2 UDs tend to penetrate 
somewhat higher into the umbral photosphere than their 
counterparts from PI. 

The elevation of the surfaces of constant optical depth 
above the upflow plumes tends to hide their magnetic 
and flow structure from s pectroscopic observations (e.g., 
iDegenhardt fc Litesl ll993a,b). It is therefore not surprising 
that observational results for the magnetic and flow prop- 



L. Bharti et al.: Properties of simulated sunspot umbral dots 



9 



3.6 
3.4 
3.2 
^3.0 
g 2.8 
2.6 
2.4 
2.2 




2600 
2400 
^ 2200 
° 2000 
ttT 1800 
§ 1600 
1400 
1200 



0.1 0.2 0.3 0.4 0.5 0.6 
UD area [Afm 2 ] 



0.1 0.2 0.3 0.4 0.5 0.6 
UD area [Afm s ] 




| -600 



2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 
mean(l)/I(DB) 



2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 
mean(I)/I(DB) 



2400 
' 2000 
1600 



400 L . i i i i i , , , i , , . , , , 

2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 
mean(l)/l(DB) 



140 
& 120 



1 00 




2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 
mean(l)/l(DB) 



Fig. 18. Relationship between properties of UDs defined 
by segmentation of continuum images at 630 nm. Mean 
quantities are averages over the area of the individual UDs. 
Red and green dots refer to datasets PI and P2, respec- 
tively; solid black (PI) and blue (P2) lines connect binned 
values for 100 points each. Positive velocity corresponds to 
upward motion, negative velocity to a downward flow. 



erties of UDs do not provide a unique picture. A reduced 
field strength in UDs has been repeatedly reported (e.g., 
IWiehr fc DegenhardtllT993l : I So cas- Navarro et al.ll200l and 
some authors also find indica tions for a decrease of 
field strength with depth (Riethmiiller et al. 2008a; 



Bharti et a l. 2009). Other studies (jSobotka fc Jurcakll2009 : 



Watanabe et al . 2009b) do not show significantly lower field 



strengths in UD s. While upwellings in the deep layers of 
UDs were found (ISocas-Navarro~ t al. 2004; RimmeTe1 l200l 
iRiethmliller et a l. 2008a; Bh arti et al. 2007a, 2009), the jet- 
like upflows predicted by the simulations in the higher lay- 
ers have not been observed so far. This may be due to the 
fact that spectroscopic studies mostly refer to large UDs, 
for which the jets are weaker or absent. 

The relations between various mean quantities for the 
full set of UDs found by segmentation of continuum images 
at 630 nm are illustrated by scatter plots in Fig. [181 The 
positive correlation between UD area and brightness seen 
in Fig. [11] (for bolometric brightness) is confirmed also for 
630 nm (top left panel). Although there is a large amount 



of scatter, UDs from P2 tend to be brighter for the same 
area, possibly reflecting higher upflow speeds (below optical 
depth unity) with a higher convective energy flux. 

This interpretation is supported by the decrease of 
the magnetic field strength with increasing UD area at 
7"630 = 0.1 (top right panel) and with increasing intensity 
at T630 = 1 (bottom left panel) for P2 as well as by the rise 
of the elevation of the level T630 = 1 with increasing UD 
intensity (bottom right panel). The downflows at tqso = 1 
intensify with increasing brightness (middle right panel). 
The upflows T630 = 1 are largely independent of bright- 
ness (middle left panel); since the level surface is higher for 
brighter UDs, this also indicates stronger underlying up- 
flows. 



5. Conclusions 

The simulation of magneto-convective energy transport in 
a strong vertical magnetic field exhibits narrow upflow 
plumes, whose properties can be compared with observa- 
tions of (central) umbral dots in sunspots. Using a multi- 
level image segmentation and tracking algorithm, we ana- 
lyzed sets of simulated UDs from two phases of the sim- 
ulation, the thermal relaxation phase (PI) and the quasi- 
stationary phase (P2). This led to the following results: 

1. Size: Histograms of UD area indicate that there is no 
'typical' size of the simulated UDs. The average area is 
0.08 Mm 2 in PI and 0.14 Mm 2 in P2. Reported values 
from recent high-resolution observations are somewhat 
lower. This is partly due to different size definitions, but 
could also indicate a lack of small, short-lived UDs in 
the simulation. 

2. Brightness: Averaged over their area, the bolometric 
brightness of the simulated UDs exceeds the surround- 
ing dark background by a factor 1.6 (PI) and 1.7 
(P2), respectively. The corresponding values for the 
continuum brightness at 630 nm are 2.6 (PI) and 2.9 
(P2). For the peak intensities, the factors can reach 
markedly higher values between about 3 (bolometric) 
and 8 (630 nm). None of the simulated UDs exceeds 
the corresponding brightness values of the quiet Sun. 
Comparison with observations is complicated by the of- 
ten unknown amount of straylight contamination. 

3. Lifetime: Average lifetimes for the simulated UDs are 
28 min (PI) and 25 min (P2). Reported lifetimes 
from observations vary significantly, but recent high- 
resolution data typically indicate shorter lifetimes, al- 
though long-lived UDs are also found. It is possible that 
the simulations, owing to limitations in spatial resolu- 
tion, miss a population of small, short-lived UDs. 

4. Correlations: Larger UDs tend to be brighter and live 
longer, although there is a significant amount of scatter. 
Similar trends have been reported from observations. 

5. Magnetic field and flows: The drastic reduction of 
the magnetic field and the strong flows in the near- 
photospheric parts of the upflow plumes are largely 
hidden from spectro-polarimetric observations. This is 
caused by the elevation of the surfaces of constant op- 
tical depth unity, which bulge upward over columns of 
hot rising plasma. Consistent with observational results, 
only a moderate field reduction and weak flow signa- 
tures are expected at optical-depth levels where relevant 
photospheric lines are formed. 



10 



L. Bharti et al.: Properties of simulated sunspot umbral dots 



In summary, the comparison of the properties of simu- 
lated and observed (central) UDs indicates that the simula- 
tions have caught the basic underlying mechanism for the 
formation of these bright structures. Differences in detail 
are not surprising, given the still unsufficient spatial reso- 
lution of both simulations and observations. Simulations 
of full sunspot s, whi ch have recently become available 
([Rempel et al.l l2009af h show that the similar magneto- 
convective processes are responsible for the formation of 
umbral dots and penumbral filaments, also clarifying the re- 
lationship between central and peripheral UDs. More com- 
prehensive parameter studies are needed to reveal the de- 
pendence of UD properties on background field strength, 
spatial resolution, and vertical extension of the computa- 
tional box. Also, the calculation of synthetic Stokes profiles 
will permit a direct comparison with spectro-polarimetric 
observations and inversions. 
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